Single-cluster dynamics for the random-cluster model 

Youjin Deng 1 , Xiaofeng Qian 2 , Henk W.J. Blote 2,3 
1 Hefei National Laboratory for Physical Sciences 
at Microscale and Department of Modern Physics, 
University of Science and Technology of China, Hefei 230027, China 
2 Lorentz Institute, Leiden University, 
P.O. Box 9506, 2300 RA Leiden, The Netherlands 
^Faculty of Applied Sciences, Delft University of Technology, 
P.O. Box 5046, 2600 GA Delft, The Netherlands 
(Dated: July 10, 2009) 

Abstract 

We formulate a single-cluster Monte Carlo algorithm for the simulation of the random-cluster 
model. This algorithm is a generalization of the Wolff single-cluster method for the (/-state Potts 
model to non- integer values q > 1. Its results for static quantities are in a satisfactory agreement 
with those of the existing Swendsen-Wang-Chayes-Machta (SWCM) algorithm, which involves 
a full cluster decomposition of random-cluster configurations. We explore the critical dynam- 
ics of this algorithm for several two-dimensional Potts and random-cluster models. For integer 
q, the single-cluster algorithm can be reduced to the Wolff algorithm, for which case we find 
that the autocorrelation functions decay almost purely exponentially, with dynamic exponents 
Zexp = 0.07 (1), 0.521 (7), and 1.007 (9) for q = 2, 3, and 4 respectively. For non-integer q, 
the dynamical behavior of the single-cluster algorithm appears to be very dissimilar to that of 
the SWCM algorithm. For large critical systems, the autocorrelation function displays a range of 
power-law behavior as a function of time. The dynamic exponents are relatively large. We provide 
an explanation for this peculiar dynamic behavior. 

PACS numbers: 05.10.-a, 05.50.+q, 64.60.-i, 75.10.Hk 



1 



I. INTRODUCTION 



[]| of the g-state Potts model ^ 



The Kasteleyn-Fortuin mapping [1] of the g-state Potts model [2j onto the random-cluster 



model provides a way to define the Swendsen-Wang [3] and related cluster Monte Carlo 
algorithms for the Potts model. These algorithms can apply nonlocal changes to the 
configuration. For systems with long-range correlations, these nonlocal methods appear to 
be very efficient in comparison with the standard Metropolis Monte Carlo method [6( which 
applies only local updates. 

The number q of Potts states appears as a continuous variable in the random-cluster 
model; the latter model can thus be seen as a generalization of the Potts model to non- 
integer values of q. There exist several ways to simulate non-integer-g random-cluster models. 
First, Sweeny applied local updates of the bond variables While the Sweeny algorithm, 
like cluster algorithms, suppresses most of the critical-slowing down, a bond update re- 
quires a nonlocal task which increases the computer time requirements. Another algorithm 
was formulated by Hu [8]. It generates percolation configurations and applies a statistical 
reweighting in order to obtain the correct averages for the random-cluster model. A cluster 
algorithm of the Swendsen-Wang type was formulated by Chayes and Machta 9| for the 
q > 1 random-cluster model. While all these algorithms, if using a random generator of 
a sufficient quality, lead to results that are only subject to statistical errors, it was found 
that the Swendsen-Wang-Chayes-Machta (SWCM) cluster algorithm (where applicable, i.e. 
q > 1) is much more efficient than the reweighting method 10]. It was also found [h| to be 
more efficient than the Sweeny method, although the latter result depends strongly on the 
sophistication of the bond update method. Detailed studies of the dynamic critical behavior 
of the SWCM and the Sweeny algorithms can also be found in Ref. 111. 

In this work, we present a single-cluster algorithm for the random-cluster model with real 
q > 1. This algorithm thus has elements in common with the Wolff m as well as with the 

CI 

SWCM cluster algorithm |9[. Since, for integer g, the Wolff method is about as efficient as 
the Swendsen-Wang algorithm, one might expect that the same holds for our single-cluster 
algorithm in comparison with the SWCM algorithm. A test of this expectation is also 
included in the present work. 

Section [III provides an explanation of the theoretical aspects of the new algorithm. We 
include a simple example of such an algorithm, and prove that the condition of detailed 
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balance is satisfied. Furthermore, we describe other variants of the algorithm that are 
applicable to models with q > 2, and describe how the algorithm can reduce to the Wolff 
algorithm for integer q. In Sec. II I II we test the validity of the algorithm, and determine its 
dynamic exponent for two-dimensional random-cluster models on the square lattice, using 
several values of q. The generic dynamical behavior appears to be very different from that 
of the Wolff algorithm. We conclude with a discussion of our findings in Sec. HVt which also 
includes an explanation of the mechanism responsible for the unusual dynamical behavior. 

II. ALGORITHM 

A. The Kasteleyn-Fortuin mapping 

We recall the mapping of the Potts model and the random-c 
bond as well as site variables, as described in Refs. and 1C 
model, see, e.g. Ref. |l2j. The Potts partition sum is 

N q 

^=[nEin ex p( X( w> 

1=1 <7; = 1 (ij) 

where the summations are on all site variables <7j, where i labels the lattice sites. The 
second product sign indicated by (ij) is on all nearest-neighbor pairs The coupling 

K is reduced, i.e., it includes a factor 1/ksT. We consider the case of ferromagnetic cou- 
plings K > 0. The Kasteleyn-Fortuin mapping of Eq. (pQ) on the random-cluster model [l| 
introduces bond variables by = or 1 between all neighbor pairs after which the site 

variables = 1, 2, • • • , q can be summed out so that only the bond variables remain as the 
degrees of freedom of the random-cluster model. Bonds = 1 (0) are considered to be 
present (absent). 

The random-cluster partition sum thus assumes the form 

1 n c 

z ° = z > = [ n e i ^ = e n ^ > ( 2 ) 

(ij) b tj =Q {b} k=l 

where u = e K — 1 is the temperature-like parameter, = ^2 hj denotes the number of 
present bonds. The number of clusters is denoted as n c . The sum on {b} is shorthand for 
the sum on all configurations of bond variables, and n b is the number of nonzero bonds in 
the k-th cluster. 



uster model on a model with 
. For a review of the Potts 
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For q > 1 one can divide the cluster weight q in two positive contributions 1 and q — 1. 
The first contribution can be associated with one of the original Potts states. To this purpose 
we introduce 'color' variables tk = or 1 for each cluster k — 1, 2, • • • , n c : 

n c 1 , r 

z »=EnE""™M(«-i)«"n p) 

{6} fc=H fe =0 

Clusters of color and 1 have weight q — l and 1 respectively. The sum on the cluster colors 
is replaced by a sum over N site-color variables £j = or 1, together with a factor d> t *|. (with 
the convention 0°=1) for each bond variable, so that each cluster contains only sites of one 
color: 

Z b = Z tb = E n^)** fl(<7 - i) 1 "*^ > ( 4 ) 

{<} {6} (y) fc=l 

where the color of the k-th cluster is denoted t s ^) where s(k) is a site in that cluster. In a 
site configuration {t} we identify 3 types of bonds (ij): 



type 


: ti 


= h 


= 


type 1 


: U 


= tj 


= 1 


type 2 


: U 


+ tj 


= 1 



Summations and products involving only one of these types of bond are specified by ap- 
pending corresponding superscripts to the pertinent summation and product signs: 

^•EEWinb-'ilW^] 

M W (ij) k=l {b} (ij) 

{b} (ij) 

where the clusters of color are labeled 1, 2, • • • , n^- The type 1 and 2 sums can now be 
executed. After rewriting the type-0 sum one obtains the partition sum expressed in site 
variables and type bond variables: 

n (0) 

z« = z»i - E E (0) DQ(9 - Q"^] ll? 1 ^ 1 + «)] • ( 6 ) 

{t} {6} fc=l (ij) 

Eq. (|SJ) specifies the probability distribution of a system of site variables ti = 0, 1 and bond 
variables bij between nearest-neighbor sites of type 0. Each term in the second sum in 



Eq. (jSJ) specifies a cluster decomposition V({b}) of the sublattice formed by sites k with 
tk = 0. Different sets of bond variables {b} may still correspond with the same cluster 
decomposition. Thus, if we replace the sum on {b} by a sum on all cluster decompositions 
of the color-0 regions, we have to insert a summation on all {&} that are consistent with T>\ 



B. The simplest form of the algorithm 

Eq. ([7]) can serve as the basis on which a single-cluster Monte Carlo algorithm can be 
constructed. This algorithm is applied as follows to a mixed configuration specified by the 
site variables £, and a cluster decomposition T> of the color-0 sites. An initial configuration 
can, for instance, be obtained from a random-cluster configuration and assigning color 1 to 
each cluster with probability 1/q. Then, a single-cluster step is executed as follows: 

1 Choose a random site i. The action taken by the algorithm depends on the color variable 



2a tj = 1, do with probability Pi = {q — 1)/? the following: form a random cluster around 
site i with bond probability p = u/(u + 1) between sites of color 1. The sites j in the 
newly formed cluster are assigned color (i.e. tj = 0) and the number of clusters 
of color is thus increased by 1. 

2b ti = 0, do with probability p2 = 1/q the following: assign color 1 to all sites of the 
cluster containing site i, and thus decrease the number of clusters of color by 1. 

C. Proof of detailed balance 

The proof of detailed balance can be formulated as follows. Consider two mixed configu- 
rations 5*i and S2, which differ only in a region C whose sites j have color tj — 1 in Si, and 
whose sites belong to a single cluster in S^, and thus have color tj = 0. According to the 
rules given in the preceding subsection, the transition probability to move from Si to S*2 is 





{V} {b}\V k=l 



U. If 
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where N c is the number of sites in region C; N is the total number of sites in the system; 
{b} stands for the n nn bond variables on the edges between nearest-neighbor sites in C; the 
combination on {b}\C indicates the sum on all configurations {b} that connect all sites in 
C into a single cluster; denotes the number of nonzero bond variables in {&}; n p is the 
number of bond variables connecting sites sites inside C with those outside C of color 1 
(i.e., the number of bonds along the boundary of C that is broken when the color of C is 
changed). The prefactor (q — l)N c /qN describes the probability that the cluster formation 
starts within C. Each of the 'present' bonds contributes a factor u/(u + 1), and each of 
the n nn — fif, 'absent' bonds a factor l/(u + 1). Also each 'broken' bond along the perimeter 
of C contributes a factor l/{u + 1). 

The rules given in the preceding subsection also specify the probability of the inverse 
move, namely from S2 to Si, as 

rd,2) = ^. (9) 

The condition of detailed balance requires that the transition probabilities T(2, 1) and T(l, 2) 
are related to the equilibrium probabilities -P(l) and -P(2) of configurations 1 and 2 respec- 
tively: 

T(2,1)/T(1,2) = P(2)/P(1). (10) 

Since the probabilities -P(l) and -P(2) are proportional to the configuration weights specified 
by Eq. ([7]), we may write 

P(2)/P(1) = W(2)/W(1), (11) 
where the weights associated with region C in Eq. (j7|) are 

W{1) = (1 + u ) n v +nm (12) 

and 

W(2) = (q-l)J2^ b - (13) 

{b}\C 

From Eqs. ([8]) and (Q, and from Eqs. (Tl2|) and (Tl3l) . we conclude that 

T(2, 1)/T(1, 2) = (q - 1)(1 + u y np ~ nnn ^2 uHb = W(2)/W(l) , (14) 

{b}\C 

which shows that the condition of detailed balance, Eq. (fit)]) , indeed is satisfied. 
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D. Other versions 



The probabilities p\ and p 2 in Sec. Ill Bl can be chosen differently, depending on the value 
of q. For 1 < q < 2 we may take px — q — 1 and p 2 — 1. For q > 2, this is not possible 
but other possibilities arise. One can generalize the algorithm by allowing more than two 
values of the color variables t{. The most obvious way is to allow n = [q] (the integer part 
of q) values with weight one, and one special value with weight q — n. Sites of the latter 
color are divided in clusters (just as before); sites of the n remaining colors are not. A 
cluster step starting from a randomly chosen site can then be specified as follows: if that 
site belongs to a cluster (thus, of the special color 0), then the cluster is erased and its sites 
are given a random color 1,2, ■• • , n with probability 1/n each. If the cluster step starts 
from a randomly chosen site of color 1,2, •■ ■ , n, then a single cluster is formed. Its sites 
receive one of the n — 1 other weight-1 colors with probability (2n — q)/[n(n — 1)] each, and 
the cluster receives the special color with probability (q — n)/n. This choice of probabilities 
satisfies detailed balance and maximizes the probability of a cluster flip. For integer q it 
reduces to the Wolff algorithm. 



E. Test of the algorithm 

We tested the single-cluster algorithm for the cases of the q = 2, 3, and 4 Potts model 
on the square lattice, by comparing its numerical results to those of the Wolff algorithm. 
We set n = q — 1 (see Sec. IIIDj) and the weight of the color-0 clusters is thus q — n = 1. 



Simulations were performed on L x L lattices with periodic boundary conditions. After 
each single cluster step, we sampled various quantities, including the densities pi of Potts 
variables in states i = 1, 2, ■ • • , q, and the single cluster size S. The single cluster size is 
counted as the total number of lattice sites in the cluster as constructed by the algorithm. 
If the number q of Potts states is an integer, the squared Potts magnetization density m 2 
can be expressed in the densities pi as 



m 2 



AED".-«) 2 = 73TB«-V#. (15) 



The sum on the right-hand side of this equation contains q terms whose expectation value 
is equal, due to the Potts symmetry. Thus, for the expectation value (m 2 ) of m 2 we may 
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write 



n 




n(q — 1) 



(16) 



with 1 < n < q. Thus it is sufficient to sample (p.; — 1/q) 2 in order to obtain (m 2 ). While 
q is taken to be an integer in this subsection, Eq. ( fT6j) still applies for general q > 1. If g 
is not an integer, n will usually be chosen as n = [q] where [q] denotes the integer part of 
q. Although, in the case n < q, Eq. (TIB]) still leads to the same expectation values as those 
obtained by averaging on the basis of a full cluster decomposition, the autocorrelations of 
m 2 may depend on the sampling method and thus be different in both cases. 

As should be expected, for Potts models with integer values of q, the Wolff and the 
present algorithm did indeed yield mutually consistent results. This is demonstrated by the 
data for (m 2 ) and (S) in Table [I] obtained by the two algorithms for the case q — 2, n — 1. 
Furthermore, since the probability to hit a cluster is equal to its relative size, it follows that 
the two expectation values (m 2 ) and (S) are equal. Our simulation results were also in a 
good agreement with this relation, as illustrated by the data shown in Tabled for the critical 
Ising model. 

Since both simulations involve the same number of samples, the statistical uncertainties, 
shown between brackets in Table HJ reflect the relative efficiency of the Wolff and the single- 
cluster algorithm. For size L — 8, the Wolff method is about 10 times as efficient as the 
present algorithm, while this difference increases to a factor of about 100 for L = 32. It thus 
appears that the two algorithms have different dynamic exponents. 

III. DYNAMIC EXPONENTS 

A. Autocorrelation functions and autocorrelation times 

Consider an observable O, whose evolution in time t' is described by the time-series 
0(t'), where each unit of t' corresponds to one step of the single-cluster algorithm. The 
autocovariance function of O is defined to be 



C (t') = (O(0)O(t')) - (O) 2 , 



(17) 



and its autocorrelation function is 



Ao(t') 



Coit') 



(18) 



C o (0) ' 
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We then normalize time t' as t = t'S/ L 2 so that the time unit of t is the average number 
of cluster steps in which each lattice site is visited once. From Ao(t) we then define the 
integrated autocorrelation time as 



1 v 

Tint.O = g + E A °(t) , (19) 
t=l 

and the exponential autocorrelation time as 



Texp '°^febg^(t)- (20) 
Finally, the exponential autocorrelation time of the system is defined as 

r cxp = sup r exPt0 , (21) 
o 

where the supremum is taken over all observables O. This autocorrelation time measures 

the decay rate of the slowest mode of the system. All observables that are not orthogonal 
to this slowest mode satisfy r exPi e> = r cxp . 



B. Integer q 

For integer q = 2,3,4, we may set n = q, in which case the color-0 clusters have zero 
weight and are thus absent, so that the single-cluster algorithm reduces to the well-known 
Wolff algorithm [4]. Such Wolff simulations were performed at criticality. The system sizes 
were chosen as powers of 2 in the range 4 < L < 4096 for q = 2, 4 < L < 2048 for q = 3, 
and 4 < L < 1024 for q = 4. Samples were taken at intervals of one single-cluster step. The 
number of samples taken for each system size is shown in Table [III 

After a fast initial decay, the autocorrelation functions for S and m 2 decay approximately 
exponentially, but with an amplitude proportional to a power of the linear size L. Except 
for the initial decay, the behavior can be described as 

A (t) oc L-^ e ~' /Tcxp(L) , (22) 

with O = S or m 2 , which implies that z- mtj o — z ex P — so- Accordingly, a data collapse is 
obtained by plotting the quantity L Ss A s versus t/r exp . This is shown in Fig. [IJ with the 
exponent of L fixed as s s = 0.37. 
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Correlations between subsequent Wolff steps are thought to arise from overlap between 
the two pertinent clusters. The average Wolff cluster size, relative with respect to the size 
L d of the system at criticality, scales with L as S oc L 2yh ~ 2d , where is the magnetic 
exponent and d — 2 is the spatial dimensionality. The probability that two subsequent 
clusters overlap may thus be crudely estimated as L id ~ iyh . The histogram of the cluster-size 
distribution is however very wide with a large-size cutoff that scales as L d ~ Vh . Since large 
clusters contribute more to the autocorrelation function than small ones, one may expect 
that the correlations at short times scale with L instead as L~ Ss with s s < 4d — Ay^. The 
results for q = 2, 3, and 4 are shown in Table IITI1 It seems that for q = 4 the Wolff algorithm 
is slightly less efficient than the Swendsen-Wang method. 

During the simulations, also the energy density E, which is defined as the nearest-neighbor 
correlation function, was sampled. The corresponding autocorrelation function is 
shown in Fig. [2] for q = 2. This figure indicates that Ae decays approximately exponentially 
as a function of t, with an amplitude that has little or no dependence on the system size. 
It thus follows that z e ^ Pt E ~ Zmt,E- The autocorrelation times r int and r exp were obtained by 
integration and least-squares fits respectively. The autocorrelation times for L > 16 were 
fitted by 

t^ e {L) = a + bL z ^ , (23) 

and similarly for r exp . The fit yields z exp ~ z intt E = 0.07 (1). This nonzero dynamic 
exponent is in agreement with the upward curvature of the data for r exp versus L on a 
logarithmic scale, shown in Fig. [3j However, we cannot exclude the possibility that the 
dynamic exponent is zero, because the data for L > 16 can also be described by T intj £(L) = 
To + lnL(ao + a\/L + 0,2/ L 2 ), which has only one more parameter than Eq. fl23l . with 
r = -1.02 (7), a = 0.76 (2), a x = 3.4 (6), and a 2 = -10 (5); this is illustrated in Fig. O 



Such behavior would mean that the Li-Sokal bound 
the two-dimensional Ising model. 



13| is sharp for the Wolff dynamics of 



C. Non-integer q 

We performed simulations of critical random-cluster systems with sizes 4 < L < 256 
for q = 1.25, 1.50, 1.75, 2.25,2.50, and 2.75, with n = [q]. Samples were taken after each 
single-cluster step, with a total number of samples of 6 x 10 7 for each L, q. The squared 
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magnetization was obtained using Eq. (|T6|) . 

The autocorrelation functions A m 2(t) and As(t) were found to display a fast decay at 
short times t ~ 0(1), then decay algebraically, and ultimately exponentially with t. Such 
a range of algebraic behavior, which extends to t ^> 1 for at large L, is absent for Wolff 
dynamics. In the case of A m 2(t), the fast initial decay at small t appears to be hardly size- 
dependent, as can be seen in Fig. HI In contrast, for As(t), the amplitude of the algebraic 
decay is found to be size dependent. 

These dynamic phenomena are very different from those for integer q, where the auto- 
correlation functions for both quantities decay almost as a pure exponential law. It seems 
that the behavior of A (t) can be described by 

A (t, L) = L~ s °t~ r °f(t/r cxp (L)) for t » 1 , (24) 

where / is a universal function. For large t, it behaves as 

f(t/r cxp (L)) ex e ^*pW with Tcxp(L ) K L ^ X p (25) 

We analyzed Ao(t,L) by attempting to collapse the data onto a single curve according 
to Eq. (jUp . The data collapses for O = m 2 and S work only approximately. This might be 
due to finite-size corrections. The results are shown in Table IIHI 

The power-law dependence of the amplitude of the exponential decay of autocorrelations 
means that, in terms of a measure of the efficiency of the algorithm, the significance of the 
dynamic exponent z exp is limited. The exponent z- mt is a better measure of the L-dependence 
of the rate of decay of correlations, because it includes the size- dependence of the amplitude 
of the decay. The unusual behavior in Fig.H]may be expected to lead to significant differences 
between z exp and z- in t. This expectation is verified by integration of Eq. (124)) which yields 
that 

r intj0 oc L Zint , z int = (1 - r a )z exp - s a ■ (26) 

Inspection of the numerical results for z exp , tq and so in Table fTTTI shows that z exp and z- int 
must have different values. For a numerical analysis of z- m %, we have, in line with Eq. ( |T9l) . 
integrated the autocorrelation functions for m 2 and S according to 

1 T 

troCT) = - + A °(t) , (27) 
t=i 
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where, presently, T assumes the meaning of a time variable. The integrated autocorrelation 
times of the q = 1.25 model for r int m 2(T) and r inti s(T) are shown in Figs.[5]and[6]respectively. 
The lines for large L are approximately straight, which reflects the algebraic decay of the 
autocorrelation functions as a function of t. As a consequence of the exponential decay at 
large t, Ti nt (T) approaches a constant. However, integration of random correlations at large 
t eventually affects the accuracy of the numerical result for Ti nt (T) so that a cutoff has to be 
applied for optimal results. For this reason, the integrated autocorrelation times for L = 256 
could not be accurately determined and were skipped from the analysis. The remaining data 
were fitted by 

r intt0 (L) = A + BL z ^o, (28) 

where A and B are unknown constants. The fits for q = 1.25 yield z intjTn 2 = 1.4 (1), and 
Zint,S — 11 (!)■ Table HVl includes the results for the z- int ,o for several other values of q. 

IV. DISCUSSION 

As stated in Sec. (TJ one might expect that the present single-cluster algorithm would have 
a dynamic exponent that is about the same as that of the SWCM algorithm |9|]. However, 
after comparing the dynamic exponents of both algorithms, we find that this expectation is 
not justified for noninteger q. The single-cluster algorithm formulated in this work represents 
a new dynamic universality class. Finding the reasons behind this curious fact should help 
us to better understand from where critical slowing down arises, and tell us something about 
how one can further develop efficient Monte Carlo algorithms in statistical physics. 

The single-cluster algorithm described above is obviously related to the Wolff [4] algorithm 
as defined for integer-g Potts models; it can reduce to the Wolff method if q is an integer. 
On the other hand, it is different in the sense that the single-cluster algorithm acts on a 
mixed configuration of site variables and random-cluster variables. 

This mixture of different types of variables is essentially the reason that the present single- 
cluster algorithm is relatively slow. In this algorithm, a number of lattice sites belongs to 
random clusters of type with weight q — [q], while the remaining sites are decorated with 
a Potts variable in one of [q] Potts states. 

As described in Sec. Ull the only process that can change a type-0 cluster back into 
an integer spin state, depends on the random selection of a site in that cluster in the 
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beginning of each cluster step. Thus, large clusters of type are short-lived and small ones 
are long-lived. It is illustrated in Fig. [7J that the single-cluster distribution for the case 
q = 2 displays a wide range of algebraic decay and an additional maximum at large cluster 
sizes of order L Vh , preceding a rapid decay at even larger sizes. The distribution shown 
in Fig. [7J represents a time average. Individual cluster decompositions deviate because of 
thermal fluctuations. The lifetime of these deviations will naturally depend on the cluster 
size. The smaller the type-0 clusters are, the longer they will persist, and this will be 
reflected in the decay of the autocorrelation functions. The pronounced maximum in Fig. [7J 
at S ~ L Vh can thus be associated with a rapid initial decay of the autocorrelations. Once 
the largest clusters of type are updated, some autocorrelations are still persisting due to 
the thermal fluctuations of the numbers of smaller clusters that remain to be updated. After 
t' single-cluster steps, the autocorrelations of the numbers of clusters with sizes S > L 2 ft' 
will be strongly reduced, while the clusters with sizes S < L 2 /t' will mostly be unaffected. 
Since the cluster-size distribution decays algebraically in a range of S, it is natural that 
autocorrelations associated with clusters that are not yet updated display a corresponding 
power-law decay in time, as long as the smallest clusters survive. After a number of steps 
of order L 2 also the clusters of size 1 will be updated. This somewhat qualitative reasoning, 
which neglects any persisting correlations after all clusters are visited, would mean that the 
longest autocorrelation time, expressed in single-cluster updates, is of order L 2 , after which 
the autocorrelations will decay exponentially. Expressed in units of t as defined in Sec. IIII A[ 
this corresponds with autocorrelations scaling as L 2yh ~ 2 at criticality. Our numerical results 
suggest that the dynamic exponent is slightly larger, namely z exp ~ 2, but the data do not 
allow a more firm statement. 

The persistence of the smallest clusters during a time of approximate order L 2 leads to a 
long time "tail" during simulations using the single-cluster method. It is this effect that we 
hold responsible for the relatively large dynamic exponent z cxp of the single-cluster method. 
However, the amplitude of the algebraic decay of the autocorrelation functions still depends 
with a factor L~ s ° on the system size L. A positive value of so therefore means that the 
critical slowing down is less severe than suggested by the value of z exp , in agreement with 
the smaller values of z- int as shown in Table HVl 

Nevertheless our findings indicate that the single-cluster algorithm, apart from displaying 
interesting dynamic behavior, is not an efficient tool to investigate the two-dimensional 
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random-cluster model. In higher-dimensional systems we have similar expectations. But 
there still seems to be a possibility that a number of single-cluster steps alternating with a 
full-cluster decomposition, which takes advantage of the fast initial decay of autocorrelations 
of the single-cluster algorithm as well as of the absence of a long-time tail in the SWCM 
cluster algorithm, will be relatively efficient in higher-dimensional systems. 
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FIG. 1: (Color online) Data collapse of the autocorrelation function of the single-cluster size, shown 
as L Sa A$ on a logarithmic scale, versus t/r exp , with s s = 0.37. These results apply to q = 2 Wolff 
dynamics. 
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FIG. 2: (Color online) Data collapse of the autocorrelation function Ae vs. t/r^^E for q = 2 
Wolff dynamics. The system sizes L are shown in the figure. The statistical uncertainties become 
appreciable at large times. 
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FIG. 3: (Color online) Semi-logarithmic plot of the integrated autocorrelation function r^t e(L) 
versus L for q = 2 Wolff dynamics. The error bars are of the same size as the data points. The 
solid (green) line is obtained from the logarithmic fit. The difference with the power-law fit would 
not be visible on this scale. The straight dashed line represents pure logarithmic behavior r oc In L, 
and serves only for the purpose of illustration. 
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FIG. 4: (Color online) Autocorrelation function A m 2 for q = 1.25 versus time t, using logarithmic 
scales. These data apply to the single-cluster simulation of the q = 1.25 random-cluster model. 
The straight line is only for the purpose of illustration, and has slope —0.25. 
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FIG. 5: (Color online) Indefinite integral r int m 2(T) of the magnetic autocorrelation function A m 2(t) 
over the time interval < t < T. These results apply to the single-cluster simulation of the q = 1.25 
random-cluster model. 




FIG. 6: (Color online) Indefinite integral Ti n t s(T) of the autocorrelation function As(t) for the 
single-cluster size over the time interval < t < T. These results apply to the single-cluster 
simulation of the q = 1.25 random-cluster model. 
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FIG. 7: (Color online) Data collapse for the single-cluster distribution P(S) as a function of the 
cluster size S for the critical q = 2 random-cluster model. The dashed line illustrates the asymptotic 
slope —2/yh = 16/15 which applies to 1 << S « L Vh . Data are shown for system sizes L = 16, 
32, 64, 128 and 256. The quantity P represents the probability that a randomly chosen site belongs 
to a cluster of size S. 

TABLE I: Simulation results for the average squared magnetization (m 2 ) and the single-cluster 
size S for the critical q = 2 random-cluster model, as obtained by the Wolff (W) and the present 
single-cluster algorithm (S) with n = q — 1 as defined in the text. The parameter L specifies the 
linear system size. The number of samples per system size is 4 x 10 6 for each simulation, and the 
number of clusters formed between subsequent samples is 2 for L < 24 and 3 for L = 32. The 
numbers between brackets show the statistical error margins in the last two decimal places. 



L 


8 12 16 20 24 32 


m 2 (W) 
m 2 (S) 


0.64693 (18) 0.58581 (18) 0.54537 (17) 0.51584 (16) 0.49305 (17) 0.45874 (14) 
0.6478 (6) 0.5861 (8) 0.5442 (9) 0.5164 (10) 0.4932 (13) 0.4610 (12) 


S (W) 
S (S) 


0.64666 (18) 0.58581 (18) 0.54544 (17) 0.51594 (16) 0.49311 (17) 0.45878 (14) 
0.6470 (6) 0.5860 (8) 0.5441 (9) 0.5165 (10) 0.4932 (13) 0.4610 (12) 
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TABLE II: Lengths of the Wolff-type simulations in Sec. IIII Bl for L > 16 and q = 2, 3, and 4, in 
units of 10 7 samples. 



L 


16 


32 


64 


128 


256 


512 


1024 


2048 4096 


Q 


= 2 


12 


12 


12 


12 


16 


16 


32 


8 8 


Q 


= 3 


4 


4 


4 


4 


8 


8 


12 


12 


q 


= 4 


8 


12 


20 


32 


48 


72 


64 





TABLE III: Single-cluster dynamics for several values of q. The exponents r s , and r m are those in 
Eq. (|24"|) for S, and m 2 , respectively; the same labeling applies to s s and s m . The values of s m are 
not significantly different from zero for noninteger values of q. For the purpose of comparison, the 
last column shows results 



HI] for uster dynamics. Furthermore, some 

data are included for integer values q = 2, 3 and 4, with the choice n = q, so that these results 
apply to the Wolff algorithm. 



q 


Ss 


r s 




I'm 


^exp 


z exp (SW) 


1.25 


0.25 (2) 


0.25 (2) 


0.00 (2) 


0.25 (1) 


2.0 (2) 


0.00 


1.50 


0.19 (2) 


0.19 (2) 


0.00 (2) 


0.19 (1) 


2.0 (2) 


0.00 


1.75 


0.14 (2) 


0.14 (2) 


0.00 (2) 


0.14 (1) 


2.0 (2) 


0.06 (1) 


2.25 


0.26 (2) 


0.15 (2) 


0.00 (2) 


0.14 (1) 


2.0 (2) 


0.24 (1) 


2.50 


0.22 (2) 


0.10 (2) 


0.00 (2) 


0.12 (1) 


2.0 (2) 


0.31 (1) 


2.75 


0.17 (2) 


0.06 (2) 


0.00 (2) 


0.10 (1) 


2.0 (2) 


0.42 (2) 


2.00 


0.37 (2) 




0.14 (2) 




0.07 (1) 


0.14 (1) 


3.00 


0.34 (2) 




0.05 (2) 




0.521 (7) 


0.49 (1) 


4.00 


0.25 (2) 




0.00 (2) 




1.007 (9) 


0.93 (2) 



20 



TABLE IV: Dynamic exponent z- mt of the single-cluster cluster algorithm. This exponent describes 
the scaling behavior of T iri t, the integrated autocorrelation function. For a negative exponent z; nt , 
the Ti n t data approach a constant when L — » oo. The values of z* nt q are calculated from Eq. (|26p 
and Table HTI| while those of ^mt.o follow from the fits using Eq. (|28j) . Some data are included for 
integer q; these results apply to the Wolff algorithm. 



q 


1.25 1.50 1.75 2.25 2.50 2.75 


2 3 4 


^int,m 2 
Z mt,m 2 


1.4 (1) 1.5 (1) 1.6 (1) 1.9 (1) 1.9 (1) 2.0 (1) 

1.5 (2) 1.6 (2) 1.7 (2) 1.7 (2) 1.8 (2) 2.0 (2) 


-0.16 (2) 0.485 (7) 1.005 (9) 


Zint,S 
Z int,S 


1.1 (1) 1.1 (1) 1.3 (1) 1.6 (1) 1.7 (1) 1.8 (1) 
1.3 (2) 1.4 (2) 1.6 (2) 1.4 (2) 1.6 (2) 1.7 (3) 


-0.4 (1) 0.16 (4) 0.72 (5) 
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